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We present a new formalism, together with efficient numerical methods, to directly calculate the 
CMB bispectrum today from a given primordial bispectrum using the full linear radiation transfer 
functions. Unlike previous analyses which have assumed simple separable ansatze for the bispec- 
trum, this work applies to a primordial bispectrum of almost arbitrary functional form, for which 
there may have been both horizon-crossing and superhorizon contributions. We employ adaptive 
methods on a hierarchical triangular grid and we establish their accuracy by direct comparison with 
an exact analytic solution, valid on large angular scales. We demonstrate that we can calculate the 
£ — , full CMB bispectrum to greater than 1% precision out to multipoles I < 1800 on reasonable compu- 

tational timescales. We plot the bispectrum for both the superhorizon ('local') and horizon-crossing 
('equilateral') asymptotic limits, illustrating its oscillatory nature which is analogous to the CMB 
OA , power spectrum. 

INTRODUCTION 

Measurements of the cosmic microwave background (CMB) radiation are now strongly constraining the viability 
of a variety of cosmological theories for the formation of large-scale structure. While the simplest 'vanilla' models of 
single- field inflation appear sufficient to explain current data [![, forthcoming experiments will test this concordance 
further. One area in which we can expect significant improvement will be in estimates of non-Gaussianity. Single- 
ts, i field inflation predicts that variations in the CMB should follow simple Gaussian statistics to high accuracy. On 
the other hand, more natural inflation scenarios which involve multiple fields and more complex potentials can, in 
principle, produce large non-Gaussian signatures that may become measurable in future experiments (see, for example, 
refs. @, HI HI, @] ) • Of course, the CMB power spectrum is not sensitive to these non-Gaussianities so we must look 
to other higher order correlators. The bispectrum is the three-point correlator of the multipoles and we note that 
statistics based on it have been shown to be optimal for testing non-Gaussianity for a variety of inflationary scenarios 
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The aim of this paper is to provide a framework for evolving the primordial bispectrum (the bispectrum at the 
■£j . end of inflation) until it leaves an observable imprint in the CMB today. Current theory dictates this can be done 
by solving a four dimensional integral over the primordial bispectrum multiplied by six highly oscillatory functions. 
In practice this is not possible to evaluate numerically. In section 2 we detail the current most popular approach 
Q which assumes a simple separable form for the bispectrum, so that the integral separates and becomes tractable. 
These methods are good for obtaining quick estimates but their applicability is limited by their ability to fit a suitable 
analytic approximation for the shape of the primordial bispectrum. In section 3, we detail a different more general 
approach. As the primordial bispectrum is the most important ingredient in the calculation, we place no restrictions 
on it other than an overall separable scale dependence. From here we reduce the four-dimensional integral to a two- 
dimensional integral over the primordial bispectrum multiplied by two one-dimensional integrals, one model dependent 
and one purely geometric. 

In section 4, we investigate analytic solutions for the bispectrum. There are two asymptotic categories for the shape 
function - "local" and "equilateral" - depending on whether the non-Gaussianity is created primarily on superhorizon 
scales or at horizon-crossing, respectively. In the large angle approximation and for the simplest form of the local 
bispectrum, there exists an exact solution which can be used to test the accuracy of our numerical methods. While 
there is no similar solution for the equilateral scenario, we argue that our methods should be more precise in this case. 
We have developed code that will accurately evolve the bispectrum from the end of inflation until today, evaluating 
its angular counterpart in the CMB. Note that our primary focus here is on the effect of primordial non-Gaussianity 
on the CMB, so it suffices to use linear radiation transfer functions for this purpose. In section 5, we detail the 
inner workings of the code and point out the key innovations required to compute the integral quickly and accurately. 
In section 6, we test against the local analytic solution and demonstrate an accuracy of less than 1% for multipole 
values less than I < 550 (using the idealized large angle transfer function). We also study the bispectrum for the 
equilateral case and find it has better convergence with an error of less than 1% for I values, I < 1300. With the 
correct radiation transfer functions used at small angles, we find that calculations of the full CMB bispectrum are 
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numerically tractable at large I and that the errors remain below 1% out to multipoles I < 1800. We show that 
growing errors for even larger I are linked directly to the resolution and the asymptotic cut-offs employed. If the 
parameters governing these are tightened, arbitrary accuracy can be achieved though at some computational cost. 
We provide plots of the CMB bispectrum for the full radiation transfer function for both the local and equilateral 
cases, showing oscillatory behaviour analogous to the angular power spectrum. Having established the accuracy and 
efficiency of these methods, we detail future directions for this work in section 7. 



BISPECTRUM 



We will quickly review the calculation for the bispectrum today from a given primordial bispectrum, this section 
loosely follows the derivation in [loj]. The temperature fluctuations in the CMB can be decomposed into spherical 
harmonics, Yi m 

— ( fi ) = ^a ;m y ;m (n). (1) 

Ira 

We can invert this expression to find the ai m 's using the spherical harmonic orthogonality condition, 



aw, 



J dn^(n)y^(n). (2) 



The temperature anisotropy can be decomposed into Legendre polynomials Pi(/J.) in k space, 

AT f d 3 k 00 

— (*)= J ^^(-ij'C^ + lJ^WAKfe^^-fl), (3) 

where 4 r (k) is the primordial gravitational-potential perturbation and A;(fe) is the radiation transfer function. We 
replace the Legendre polynomials using the spherical harmonic addition theorem, 

4-7T ^ 

P t (k . fi) = — £ Y lm (k)YC m (h). (4) 

m=—l 

Making these substitutions, and using the orthogonality of the spherical harmonics, we have, 

d 3 k 



d 3 n £(-<)'' (21' + l)tt(fc)A,, (k)P v 0*)Y? m {n) 



(2tt)3 
d 3 k 



J ^E(-0''*(*)^(*)X;i7 m <(k) / <PnY Vm ,(n)Y? m (n) 



= M-*) 1 I j^vWMWmfc) (5) 



d 3 k 

W) 

The bispectrum is defined as the three-point correlator of the a; m 's, 



hhh — \ a iimi "i2m2 a (3"i3 / v D / 

_ (4 . ) .(_ i) , + .. + ., / M«WM«*)) 

A,, (fa ) A,, (fa)A,, (*,)y,; mi (feW™, (fc)l',:,,,, (is). (7) 
where fci = |ki|, fc 2 = |k 2 | and fc 3 = |k 3 |. The primordial bispectrum is defined as 

(*Cki)¥(k 2 )tf (ka)) = (2n) 3 F(k 1 ,k 2 ,k 3 )5(k 1 + k 2 + k 3 ) , (8) 
where F(k\, fc 2 , k 3 ) is the primordial bispectrum shape. We use the integral form of the delta function, 
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and we replace the exponential with spherical harmonics using the Rayleigh expansion, 

l m 

where x = |x| and the unit vector x = x/x. Substituting, then using the orthogonality of the spherical harmonics 
and the Gaunt integral, 



G hhh .hoy Y Y - , / {2h + 1){ ~ 2h + 1)(2 * 3 + l ) (hhh\( h h h , 
y mi m 2 rn 3 - J anr hmi r l2m2 Y hm3 - y ^ \Q Q Q J \ mi m 2 m 3 I ' ' 

' h h h \ 

where ( I is the Wigner 3j-symbol. We find the bispectrum becomes 



B u2 2m3 = (-) G l ^[±m s I ^1*2*3 (zfeifc 2 fc 3 ) 2 ^(^^ 



mi mi nT,3 

'1,1 w\ ~~ 1 ~ ) t',,, ,„_,,„ . 

It is worth noting that the Wigner 3j-symbol places some restrictions on the I values. First, no individual I can be 
greater than the sum of the remaining two (the triangle condition) and, secondly, that the sum of all three must be 
even. As we assume statistical isotropy it is more common to work with the angle averaged bispectrum, 

V ( h h h \ B™ 1 ? 2 " 13 - (13) 
t-^ \ mi m 2 m 3 I llt2t; > y ' 
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Observing that, 



y [ h h h V = i, (w) 

^-^ \ Till 1Tl 2 m 3 / 
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we have, 



2\ 3 /(2/ 1 + l)(2I 2 + l)(2/ 3 + l) (h h h 



4tt V 

dxdkidk 2 dk 3 (xkik 2 k 3 ) 2 F(ki, k 2 , k 3 )A h (ki)A h (k 2 )A h (k 3 )j h (hx)ji 2 {k 2 x)ji 3 {k 3 x). 



(15) 



Also noticing that the Gaunt integral represents a purely geometrical factor which is independent of the primordial 
bispectrum we define, 

hhh ~ ym 1 m 2 m 3 hhl3 ■ 

where &z 1 z 2 ; 3 is the reduced bispectrum, 

2 x3 



h hhh = - dxdhdk 2 dk 3 (xkik 2 k 3 ) F(ki,k 2 ,k 3 )Ai 1 (ki)A[ 2 (k 2 )A h (k 3 )ji 1 (kix)ji 2 (k 2 x)ji 3 (k 3 x). (17) 



This is the quantity that we will be calculating. 

To proceed in general we must solve a four-dimensional integral of the shape function multiplied by six highly 
oscillatory functions, A/ and ji, which is not feasible in practice. Current approaches overcome this by assuming 
a simple separable form for the shape function. The most common is that developed in ref. Q where the primary 
assumption is that the primordial gravitational potential perturbation can be approximated by, 

ttfc) « * L (x) + f NL - , (18) 

with $i(a;) the linear Gaussian part and /jvl & constant. Following this through gives a shape function of 

F(k u k 2l k 3 ) = 2f NL (P*(fcx)P*(fc 2 ) + P*(fc 2 )P*(fc 3 ) + P*(fc 3 )P*(fci)) . (19) 
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Substituting this into eg. (|17|) we can now separate the integral into, 

J x 2 dxbf i {x)bf 2 {x)bf 3 L {x) + perms (20) 

with 

bf(x) = ~ J k 2 dkP^(k)&i{k)ji{kx) (21) 
bf L (x) = -f NL /VdfcA,(fc)j,(fcr), (22) 



that is, it becomes products of one-dimensional integrals which are manageable numerically. 

Searches for the bispectrum in the CMB have consequently focused on establishing limits on f^L via purpose-built 
statistics rather than calculation of the bispectrum directly. The WMAP team have constrained non-Gaussianity [l[ 
by narrowing their analysis to the local ansatz (fl~8|) . estimating the limits to be —54 < Jml < 114. (This analysis 
has been repeated in ref. [8] with an improved estimator to obtain the limits — 36 < /ml < 100.) Now primordial 
non-Gaussianity can be broadly classified into either of two asymptotic regimes, "local" or "equilateral" [ll|. The 
local case describes models in which the non-Gaussianity is produced outside the Hubble radius ('supcrhorizon') with 
examples including the curvaton model 0] and multiple field inflation models The equilateral case describes the 
case in which the primary contribution to non-Gaussianity arises just as the perturbations leave the Hubble radius 
('horizon-crossing'), with examples including DBI model [H and ghost inflation Following the 'local' WMAP 
approach described above, there has also been a parallel analysis for the equilateral case [20j|. The key element shared 
by this approach is an analytic approximation for the shape function that allows the separation of the bispectrum 
integral into the product of three one-dimensional integrals over ki inside an integral over x as in (|20p . They then 
proceed by generating a statistic sensitive to their ansatz to establish limits on a generalisation of the non-Gaussianity 
parameter Jnl- They find considerably weaker constraints on /jvx m the equilateral case, —256 < Jnl < 332 0], 
but argue that this represents approximately the same "level of non-Gaussianity" . Other recent work in this area 
includes ref. [l3| where, again using the separable local ansatz (|2"0|) . the particular scale dependence of fNL generated 
by post-inflationary gravitational evolution is incorporated. We note also for the local assumption (fl8|) in ref. [l4j ]. 
map-making methods are developed to create individual realisations incorporating the bispectrum and higher order 
correlators. 

In this paper we have adopted a different more general approach which does not rely on separable ansatze giving rise 
to the triple decomposition (|20p . Rather than generating approximations to the bispectrum to allow for easy analytic 
calculation, instead we have written code that will directly solve the integral for any shape function assuming only 
the separation of an overall scale dependence. We anticipate that this analysis will encompass all viable inflationary 
models. 



NEW PARAMETRISATION 

By looking at the expression for the reduced bispectrum (|17|) . we can rearrange the terms as follows: 

dk 1 dk2dk 3 (k 1 k 2 k 3 ) 2 F(k!, k 2 , k 3 )Ai 1 (k 1 )A h (k 2 )Ai 3 (k 3 ) I / x 2 dxji 1 (k 1 x)j h (k 2 x)ji 3 (k 3 x) ) . (23) 



2 x3 



The integral in the brackets arises as a direct result of the delta function and analytic solutions have been found, 
[TBI [lil 17 1 . These solutions demonstrate that the three lengths k\,k 2l k 3 must be such as to be able to form a 
triangle, as expected from the delta function, which restricts us to the region of fc-space highlighted in figure Q] The 
analytic solutions also state that the value on the boundary is half that of a point immediately interior. It is natural 
to reparametrise this region into triangular slices [3] , 

k x = ka = k (1-/3) 

k 2 = kb = ifc (1 + a + (3) 

k 3 = kc = h(l-a + P), (24) 
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where a, 6, c are shorthand for the corresponding expressions involving a, (3. They are not independent of each other 
as a + b + c = 2. Here, the overall wavenumber scale is given by the semi-perimeter, k = ^(fci + &2 + ^3)- The surface 
k = const defines a plane with normal (1, 1, 1) at a distance -^k from the origin and a, f3 parametrise our position on 
that plane. The new parameters have the following domains < k < 00, < j3 < 1, and —(1 — (3) < a < 1 — (3. This 
parametrisation has volume element dk\dk2dk^ = k 2 dkdad/3 (representing a minor improvement over the original 
triangular parametrisation given in [Hj]). Making these substitutions we have 

-\ J k 2 dkdad(3(abc) 2 k 6 F(ak,bk,ck)A h (ak)A l2 (bk)A h (ck) (J x 2 dxj h (akx)ji 2 (bkx)j h (ckx)\ . (25) 



and, if we rescale x to absorb k, we have 

dadf3{abc) 2 (^J ^-k 6 F(ak,bk,ck)A h (ak)A h (bk)A h (ck) ) ( / x 1 dx} u (>/.r!,/;,( h.r !,/, ; (rr ) j . (2(i) 



2 x3 



Given strong observational limits on the scalar tilt we expect the shape function to exhibit near scale-invariance 
[2> EH j so we make the substitution 

k n 

Fik^k^k^tt—Fia,®, (27) 

where n is the bispectrum tilt which is expected to be small. This is the only restriction we place on the shape 
function, up until this point the discussion has been entirely general. In fact this is slightly over-restrictive and we 
should note that the method continues to work if we replace k n with a more general overall scale dependence f(k). 
Hence, we can rewrite ([23)) as 



2 x3 



dadf3F bI {a,f3)I I (a,/3)I G {a,f3), (28) 



where, 



/dk 
A h (ak)A h (bk)A h (ck)k n — (29) 

I G {a,f3) = J ji 1 (ax)j h (bx)ji 3 (cx)x 2 dx (30) 
F SI (a,f3) = {abcf F(a,/3), (31) 

where we have defined the "scale independent" shape function, F SI . (A similar concept was discussed in ref. but 
can be contrasted with the definition (131|) in our more symmetric triangular parametrisation.) 

The end result is that we have broken up the four-dimensional integral for the reduced bispectrum (fT7j) into a 
two-dimensional integral (|28p . restricted to an equilateral triangle, over the scale- independent shape multiplied by 
two one-dimensional integrals (|2"9")) - (|3TH) . The first I T {a : j3) is over the transfer functions and the second is purely 
geometric I (a,0). These two one-dimensional integrands consist entirely of highly oscillatory functions. However, 
as we shall discuss later, since they are only one-dimensional it is feasible to evaluate them accurately numerically, 
freeing us from the severe difficulties previously encountered. 

The two basic asymptotic limits (local and equilateral) for the scale-independent shape can be approximated by 
the analytic expressions [2fj| : 

FfJUaAc) = " 3 + f b + ^ (32) 

^equilateral( a ^T C ) = a~b~C ' 

We have plotted these two bispectrum shapes in figure [2] While much emphasis has been placed on these two limits, 
we note in passing the caveat that even single-field inflation yields a more complicated non-Gaussian shape [l9| . 

We will take this opportunity to define a general momentum-dependent /jvl that will work in both the local and 
the equilateral cases (as in |12l ]) 

9f (k k k\ = F(fci,fc 2 ,fc 3 ) 

/JVM 2 ' 3j P*(fc 1 )P*(fc 2 )+P*(fc 2 )P*(fc 3 )+P*(fc 3 )P*(fci) ' K ' 
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Figure 1: Region of integration for the bispectrum calculation. The red lines are ki = fej, &3 = 0; fe = ^3, fei = 0; kg = 
ki, ki = and the region of integration is in yellow. This area can be parametrised into slices represented by the green triangle 
and the distance -y=k of the centre of the triangle from the origin. 




Figure 2: F SI plotted on the a/3-triangle for the 'local' superhorizon shape (left) and the 'equilateral' horizon-crossing shape 
(right). The equilateral case has been scaled to that the centres of both plots are at the same height. 

This agrees with the definitions given in both [3] and [23], given their respective assumptions, but it can now be 
applied across the full ki parameter space. Further, if we make the same parameter substitutions as those used in 
obtaining Q28p. we see that Jnl bas an overall scale dependence given by n^L = n + 2(1 — n s ) where n s is the usual 
scalar spectral index of the power spectrum. 



ANALYTIC RESULTS 



As well as the analytic approximations for the shape function in the local and the equilateral limits (|32p . on large 
angles (I <C 200) the transfer function can be approximated by, 

A l (k) = ~j l (A V k), (35) 

where A77 is the conformal time period between today ?y an d the surface of last scattering ?/dcc- Substituting these 
into the bispectrum, the whole integrand can be expressed in closed form and so we can look for analytic solutions. 
In the local case, we make the replacement 

^ I (*i>fa.*») = iS- + T^ + inr ( 36 ) 

K2K3 k 3 ki k\k2 
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and the bispectrum integral separates into three integrals over the fej's inside an integral over x, 

J-) 3 J x 2 dx kldk l3h (A V h)j h (xk!)) (| (A^ 2 )j; 2 (xfe)) (| ^ji, (A V k 3 )j h (xk 3 )^j , (37) 

plus permutations. The k\ integral evaluates to j^i8(x — Arj) and integrating out the delta function leaves, 

.„7)(/**<H(/** ( H- (38) 

These two integrals evaluate to (21(1 + so with permutations we have the analytic result 

n\f i . i . i ^ (39) 



\27iry \h(h+l)l 2 (h + l) Uh + lMh + l) l 3 (h + l)h(h + I) 



This exact solution (in the large angle approximation) provides an excellent reference for checking the accuracy of our 
numerical methods in later sections. 

For the equilateral case, we make the replacement 

tpsifu u n 1 (-fci + k 2 + fc 3 )(fci ~k 2 + fc 3 )(fci + fc 2 - fcg) f , n , 
t (fci,fc2,K3) = r r-TT v 4U ) 

= - ( — 2 — — — — h h all possible permutations ) . (41) 

8 V k 2 J 

When we substitute this into the full integral it again separates into three integrals over the ki's inside an integral 
over x. The middle term is the same as the local case with the exact solution given above. For the remaining two 
terms we are left to solve integrals like, 

dkk n j,{k)Mxk) = 2^/ dkk n - l J l+h {k)J l+h {xk), (42) 

for n = 1, 0, — 1. The case n = has a nice solution [23| (p405), 

ttx-V+V n x l 

x>l : x<l => . (43) 

2 21 + 1 2 21 + 1 V ' 

Unfortunately, however, both n — 1,-1 produce complicated expressions involving hypergeometric functions [23| 
( P 401), 

->1 - ^^,-(-) 2 ^(l, + 1 , + |J,) (45) 

/ir(0 , _ / i , , , 3 , 



n = -l : x<l =► 4 p^ + 3^ aJi (^--, /;* + -; 2^) (46) 

so apparently there is no analogous route to a simple analytic result. 

We have also pursued an alternative approach by deriving series solutions for the integrals required for the bispec- 
trum. These are presented in the Appendix, along with a minor correction to the literature, but their utility is yet to 
be established since they exhibit slow numerical convergence in key limits. 



NUMERICAL SOLUTIONS 



The main problem with a numerical evaluation of the bispectrum is the intensive sampling required to deal with 
the extreme oscillatory nature of the integrand. Recall the progress made so far with a 4D integral over six highly 
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Figure 3: A plot of the spherical Bessel function I — 200. Note how the function essentially vanishes up to I ~ 180. 

oscillatory functions broken down into two ID integrals (|29|) and ([30]) over three highly oscillatory functions each, 
inside a 2D integral (|28j) over a fixed triangular domain. This is a significant improvement because it has reduced the 
overall dimensionality of the problem and limited the worst oscillatory behaviour to one dimension. However, we still 
have to work hard to get the integrals to converge fast enough for the calculation to be feasible. At each point in the 
a/3-triangle we must solve the two ID integrals, I T (a, j3) and I (a, (3) with each integrand composed of either Bessel 
or transfer functions which are non-trivial to calculate. Fortunately, the only part of the integrand that changes as 
we move about the a/3-triangle (i.e. at fixed l\, l 2 , h) are the values of a, b, c (see ([M])). These simply rescale the 
three functions relative to each other so we do not need to re-evaluate them at each point. Instead we can calculate 
the functions once beforehand, retaining them for future reference. At each point in the a/3-triangle, we then rescale 
and multiply them together to form the integrand, saving lengthy calculation. 

The preliminary calculations of the Bessel and transfer functions are performed for a closely spaced range of l- 
values with the results stored in a table, which can then be read in prior to every bispectrum calculation. The transfer 
functions for all Vs were calculated using the line-of-sight code CMB2000 (or else an equivalent such as CMBFAST 
or CAMB), along with the associated spherical Bessel function at the same /, and stored in separate tables. Next we 
take the tabulated functions and interpolate using a cubic spline to create a dense grid of values. Bessel and transfer 
functions for large I start with a long flat section which is approximately zero before the first peak and hence does 
not contribute to the integral (see figure [3]) . This section is set identically to zero with a view to skipping it in the 
later integration; we define this as the part for which the function magnitude remains 10 -10 times below its maximum 
value. 

Now consider the two-dimensional integral (|28p over the a/3-equilateral triangle at fixed multipoles l\, h- Having 
produced the two tables of Bessel and transfer functions we access the relevant rows for the given multipole set. At 
each point in the a/3-triangle we then calculate a, b, c and stretch the data appropriately using linear interpolation to 
calculate the new points. We skip the first section where the values are zero and then calculate the two onc-dimcnsional 
integrals using the trapezoidal rule. This stretching and linear interpolation will tend to magnify any errors in our 
original data, so this is mitigated by using a very large number of points (O(10 6 )), even though the integrands use 
a much smaller selection (say every 20th point). This balances accuracy, using dense initial sampling, with speed, 
using sparser sampling for the integration. Next we recall that for the geometric integral the value on the boundary 
is half that of the interior, causing slow convergence of the two-dimensional integral near the discontinuous edges. 
We avoid this problem by doubling the result of the geometric integration for points actually on the boundary. The 
important further check that needs to be made is when to terminate the integration. There is a hard boundary when 
the tabulated data runs out (set by an overall range), but most integrals converge well before this limit. Convergence 
is determined by comparing the contribution of the last x sections (typically 0(10,000)) to the value of the total 
integral. If it remains unchanged up to a certain percentage tolerance (typically « 0.0001%) then the integration is 
stopped. However, given the oscillatory nature of the integrand it is important to cover a large region over many 
periods in order to correctly determine this threshold. This is most important in the corners of the triangle where the 
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Figure 4: The equilateral triangle we need to integrate over to obtain the bispectrum for h = I2 = h = 20. Note the small scale 
structure to the right of the peak. The first few oscillations make a significant contribution and must be sampled intensively 
for an accurate result. The rest of the oscillations, however, do not because of cancellations. Hence any adaptive algorithm 
must be able to determine the importance of any structure it finds. 




Figure 5: Refinement method. We start with the black cell defined by the three black points at the corners. The cell is then 
divided into four by calculating the three red points and the change in area is then j ^ black — ^2red\. 

data is stretched, along with the period of oscillation, by the largest amount. As we shall note in the next section, 
these truncation errors in the ID integrals can be significant at large / unless handled carefully. 

When attempting to perform the 2D integral (|28|) over the triangle we have to balance two key considerations. 
Every point used entails the calculation of both one-dimensional integrals so we want to minimise their number to 
reduce calculation time. Unfortunately, while the majority of the highly oscillatory behaviour is confined to the 
one-dimensional integrals there is still some oscillatory behaviour which carries over to the a/3-triangle (see figure [4]). 
As the a/?-triangle has small scale structure we need to use a fine mesh to gain accuracy. To achieve this we have 
developed an adaptive method that first uses a sparse triangular grid to estimate the integral using a 2D-trapezoidal 
rule. This then steps through each triangular cell and divides it into four, calculating the change in the local integrand 
contribution (see figure [5]). If this change is below a tunable accuracy parameter then it uses the original result, but 
if it is above then it continues the refinement process by dividing each of these triangles into four new triangles and 
estimating their contribution. This will be repeated until the change in the local result is below an acceptable limit 
(or the recursion hits a predefined tunable recursion depth). The net effect is that we only have dense grids in areas 
where the most significant contributions to the integral are made, with the remaining areas only sampled sparsely. 
For the h = h = h = 20 triangle integration, we achieve convergence for a recursion depth of five after which there 
is little improvement. One must be careful, however, to set the original sampling grid to broadly capture most of the 
structure or else the refinement algorithm may cease prematurely (see figure [6]). 

We note that these calculations naturally coarse-grain the computational work either through the sampling of ID 
integrals on the 2D triangular grid or, at a higher level, simply by evaluating the bispectrum at different multipolc 
values. The problem is well-suited to parallelisation on a large supercomputer or cluster and this has been achieved 
with the present code using an MPI implementation which dramatically reduces calculation timescales. 
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RESULTS 



As we previously noted, bispectrum multipole values l\, 1%, I3 are restricted by the same triangle condition as the 
wavenumbers ki, ki, k%. This allows us to propose a similar parametrisation to that we used for the fcj's in (|24|) . that 
is, 

h = §i(l-7) 
h = |z (1 + 5 + 7) 

h = ^(l-5 + 7 ). (48) 

We will use these parameters to graph the bispectrum in two ways, (i) the equal multipole case l\ = 1% = I3, i.e. 
the distance from the origin plotting against I (holding (5 = and 7 = i) and (ii) transverse triangular slices with 
I = const. 

In order to characterise the main sources of systematic error, we begin by comparing to the exact solution in the 
local case ([3U1) using the large angle transfer functions (|55|) . This is a stringent test because the fall-off in the radiation 
transfer functions is much faster at large I in the realistic case; the errors in the test case are more pronounced and 
help delineate both resolution and truncation errors. We have calculated the equal momentum bispectrum from I = 2 
to / = 1800 on a fairly sparse initial triangular grid of sidelength N = 200 and a maximum recursive depth of only 2, 
and we find that the numerical result has an error less than 1% for I < 550 and a maximum error of 9% for I — 1800. 
We also repeated the same calculations with a much finer grid with iV = 1600 to find comparable results which are 
shown in figure [7l demonstrating little improvement from the increased resolution. The calculations were performed 
on the Cosmos supercomputer and took between 16 minutes and 4 hours 55 mins cpu time with a median of about 80 
mins for the adaptive method. The average time to calculate each bispectrum point bm on the finer grid was greater 
than 11 hours. 

The largest source of error in the calculation is a tendency to underestimate the integral due to premature truncation 
of the one-dimensional integrals for large /. There are two reasons for this. First, in the case of I = 1800 the Bessel 
function is approximately zero until about 1700. Truncating at a value of 10000 we only need to rescale by a < 0.17 
for the non-zero section of the integrand to be shifted off the interval, making the whole integral vanish numerically. 
However, this only affects the outside of the a/3-triangle so it only generates a small error. The second reason is that 
the intelligent truncation routine that decides if the one-dimensional integral has converged can be too aggressive. 
This affects all points and is the primary source of numerical error. If we double the integration cut-off and eliminate 
the intelligent truncation the error in the one dimensional integrals is drastically reduced (see figure [8]) . The cost is 
a fivefold increase in calculation times for an adaptive N = 400 grid. This is mainly a problem for the geometric 
integral as this only decays as l/x requiring late truncation for accuracy. Here, the large-angle transfer integrand 
decays as l/x 4 so it converges more quickly. If we calculate the test case again with an N = 1600 grid, doubling the 
ID region of integration, and turning off the truncation routine, then the error is reduced comfortably within 1%. 
Unfortunately the integrals then take just under 80 hours to complete. 

Having confirmed the accuracy of the method, in principle, the calculations were then repeated with the full 
radiation transfer function. Here, with an initial N — 200 grid and a recursion depth of 2, the calculations completed 
on Cosmos in a median time of 31 mins. The results for the equal momentum bispectrum are shown in figure Its 
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Figure 6: Successive refinement of the grid for integration of the h = fa = fa = 20 triangle, the left is a recursion depth of 
one proceeding to five on the far right. After a recursion depth of five there is little change in the grid or result, however the 
error for the integration is still 4.45%. If we look at the diagram we see a blank triangle about halfway between the corner and 
the centre. This region contains important structure but it has been missed as the initial sampling grid was too sparse. By 
doubling the initial grid the error is reduced to 0.8%. 
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behaviour is as expected with the two main peaks appearing at I w 200 and I ~ 500, mirroring those observed in 
the CMB power spectrum. We have also calculated the bispectrum using a dense TV = 1600 grid with the difference 
also plotted in figure O confirming the absence of a significant resolution error. On the other hand, we might be 
concerned about the important truncation error we observed in the large-angle test case discussed above. To check 
this we ran long calculations for the N = 1600 grid but without intelligent truncation and while also doubling the 
cut-off in the geometric one-dimensional integral. The results gave a maximum correction of only 0.004%, indicating 
that truncation error is also insignificant in the realistic case. The explanation for this improvement in accuracy lies in 
the much more rapid exponential fall-off at large I of the true radiation transfer function, when compared to the large 
angle solution (135|) . Since the transfer function integral multiplies the geometric integral in the final 2D integrand 
(|28[). it also suppresses any errors due to the slow convergence of the latter. The comparison in figure [9] indicates 
that even the N = 200 adaptive grid achieves a bispectrum accuracy better than 1% across the full multipole range 
I < 1800. 

In the equilateral case we no longer have an analytic solution with which to compare. However, both I T (a,f3) and 
I (a, f3) are strongly peaked in the centre of the a/3-triangle (as we can see from figure |4j. Since both bispectrum 
shapes are flat in the centre, we expect the equilateral case to exhibit approximately the same scaling with I as the 
local case for bin. To test accuracy we can again take the large angle approximation, dividing by the local analytic 
solution to normalise. As anticipated, we find that numerical accuracy is improved over the local test case, because 
the more oscillatory contributions near the boundaries are further suppressed by the primordial shape. Using an 
adaptive N = 200 grid (recursive depth 2), we find that this test case has an error of less than 1% for I < 1300 and 
maximum error of 1.7% for / = 1800 when compared to the calculation on a N = 1600 grid without the intelligent 
truncation and doubling the cut-off (see figure fTO]). 

The calculation for the equilateral case with the full radiation transfer function gives much the same results as for 
the local case (see figure ITTj) . Here we again divide by the analytic solution scaling the result by the factor 24 so 
that the centres of the two shape functions are both the same height. The main difference between the two is that 
the peaks are slightly smaller in the equilateral case for low I (figure ITTj) . For large I, the results are almost identical 
because the larger the Z-values the more peaked I T and / become when plotted on the a/3-triangle. For large I they 
act like a delta function picking out only the shape function value where k\, ki-, k% have the same ratio as l±, I2, h- 
So for large the equal I bispectrum is essentially only proportional to the height in the centre of the shape function 
rather than its shape. These bispectrum calculations on Cosmos took an average of approximately 30 mins for each 

It is clear from figures [9] and [TTJ that it would be difficult to distinguish between the local and equilateral cases 
on the basis of the equal momentum bispectrum. Instead, to achieve this we must look at the transverse triangular 
slices. For both cases three slices were made at 3/ = 850, 31 — 1850 and at 31 — 3650 with points being calculated 
at li's that are multiples of 50. For the slowly converging test case with a local shape function and the large angle 
approximation, we find that the error in the calculation for 31 = 850 is less than 0.6%, for 31 = 1850 less than 3%, 

Test / Analytic 




Figure 7: Plot of the bispectrum for the local case in the large angle approximation divided by the analytic result. The blue 
line is with the adaptive code starting with a grid of 200 and allowing two recursions, the red is for a grid of 1600. For the fine 
1600 grid the oscillations disappear and the line is smooth indicating that this eliminates error from the grid sampling. Note 
the adaptive method remains within a percent of this line. Both the results are within a percent up to I « 500 then we notice 
a downward trend due to the premature truncation of the one dimensional geometric integral 
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Figure 8: Sources of error from the one dimensional integrals. The four lines are all for the same grid of 400 on the triangle. 
The red line is the calculation for the standard conditions, the blue is when we double the region of integration, the yellow is 
with the intelligent truncation turned off and the green is when we do both. 
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Figure 9: Plot of the bispectrum in the local case using the full radiation transfer function over the result for the large angle 
approximation. Note the three main peaks where we would expect from the plot of the power spectrum. The red line is the 
difference from using the dense grid. At 650 it reaches a maximum of only 0.045 which justifies our use of the much faster 
adaptive method. 



and for 31 = 3650 less than 5%. However the error is always negative and its range is small in all cases, indicating 
that it is primarily systematic and due to the premature truncation issue discussed above. Most importantly we note 
that the error appears to be a function only of the sum, and not the specific combination, of the ij's so there are no 
new issues in dealing with the slices. 

For the local case with the full transfer function with an adaptive TV = 200 grid and recursion depth 2 the typical 
calculation time was 44 mins. In the triangular bispectrum plots (see figure [T2"]) , we have l\ constant along the lines 
parallel to the left bottom edge, where it is at a maximum (400, 900 and 1800 respectively), reducing linearly to a 
minimum in the top right corner (50 in all three cases). The dependence on I2 and ^3 is similar with respect to the 
remaining two edges. As we used the analytic version of the local shape function, which has no initial structure, the 
diagrams are easily understood from figure [H For 3^ = 850, each of the edges corresponds to one li = 400. As 400 is 
close to zero in the plot along the equal I direction, the edges are close to zero in the slice. On the other hand, near 
the centre is the triple (250, 300, 300) all of which are large and negative so they produce a broad central trough. In 
the plot of 31 = 1850 we have a three peaks at the triple (450, 700, 700) (plus two permutations) as 700 is large and 
negative and 450 is positive. We also have troughs for triples (250, 800, 800) (near the corner) and (500, 500, 850) 
(near the middle of the edge) as both 250 and 500 are close to maximums and give strong contributions. The plot 
of 3^ = 3650 has many peaks and troughs arising from various combinations of the k's but they are an order of 
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Figure 10: Plot of the bispectrum for the equilateral case in the large angle approximation over the analytic result. The result 
has been scaled so that the result for I = 20 is 1 as this is the most stable part of the plot. As we have no analytic solution 
to compare to in the equilateral case we have calculated the bispectrum on a dense N = 1600 grid without the intelligent 
truncation and doubling the cut-off. By comparing the two plots we see the error in within 1% up until I = 1300, a significant 
improvement over the previous case. 
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Figure 11: Plot of the bispectrum for the both the local case and the equilateral case using the full radiation transfer function 
over the analytic result. The equilateral case has been multiplied by 24, the ratio between the height of the centre point of the 
local case to the equilateral case. We again see the main peaks where we would expect from the power spectrum calculation 
but the local case is smaller for low I. For large I the two one-dimensional integrals act like delta functions picking out the 
point on the shape functions where ki, &2, k^ have the same ratio as h, h, h- As we have scaled the centre points to be the 
same height the bispectrum curves merge for large I. 

magnitude smaller than the previous results. This is because all triples have at least two U > 900 so the contribution 
from the remaining li is suppressed. We can see a 3D comparison in (see figure 1171) 

For the equilateral case with an initial N — 200 grid, the bispectrum calculations were completed on Cosmos again 
in about 44 mins and are also plotted in figure [12] (right). We see the same basic pattern of peaks and troughs as 
in the local case. The main difference is that the heights of these features become strongly suppressed towards the 
edges of the triangle. This is a direct reflection of the difference in the two shape function shapes. The local shape 
(superhorizon case) magnifies these features near the edge with strong correlations between very different multipoles, 
whereas for the equilateral case (horizon-crossing) they are suppressed. 

Having proved the method for the two standard analytic cases we can proceed to look at cases where the shape 
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Local Equilateral 




Figure 12: Plots of the bispectrum for the local case (on the left) and for the equilateral case (on the right) for / = constant 
slices using the full radiation transfer function. The top row is 31 — 850, the second 31 = 1850, and the third 31 = 3650. The 
bispectrum has been divided by the analytic result so its features are clearly visible. Lines of constant U are those parallel to 
their respective edge of the triangle. The extremums of the graphs always appear when the U triples contain Us which represent 
extremums of the plot along the I direction. Note that the outside of the graphs are suppressed in the equilateral case compared 
to the local case due to the differing shape functions. 



function cannot be separated. A good example is the DBI model [5j. The shape function for this model is identical 
to than in 21[ where higher derivative operators have been added to the Lagrangian. This shape function has the 
explicit form of, 
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Figure 14: F plotted on the a/3-triangle for the 'Higher Derivative' model. 



F SI (k!, h2, k 3 ) = fcifcafc (fcl+fc2 + fc) 2 E k * + E( 2fc i " 3fc ^) + E - ^ k P0 ■ ( 49 ) 

This shape function has been plotted in figure UM 

Due to the factor (fci + k^ + k ^) 2 this shape cannot be separated as previously, although some progress has been 
made by using an integral form 22]. Subsequently we have always been forced to calculate the bispectrum for this 
model by approximating with the equilateral shape. With this new approach we can now we can calculate it in full 
as the non-separability no-longer poses a restriction. 

The calculations run in similar times as for the equilateral case. The equal 1 bispectrum gives almost identical 
results as the equilateral case as can be seen in (see figure [15]) so it may seem that the previous approach is justified. 
When we plot the transverse slices however, the differences become apparent (see figure [T6f . From these plots we can 
see that the variation between the two cases is at the 10% level. A 3d comparison can be seen in (see figure ??). 

As a final comment, we note that we can exploit the properties of the bispectrum discussed above to speed up its 
evaluation further. Very substantial improvements could be achieved by tabulating the product of our ID integrations 
over the Bessel and transfer functions (j2"91 13"0)) . rather than keeping tables of the functions themselves. This would be 
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Figure 15: Plot of the bispectrum for the both the higher derivative case and the equilateral case using the full radiation 
transfer function over the analytic result. As the two shape function shapes are similar in the region surrounding the centre of 
the triangle the results are almost identical. To see the difference between the two shapes we must look to the transverse slices. 

relevant, for example, if the background cosmology was fixed (i.e. the transfer functions are not modified), perhaps 
for a likelihood analysis marginalising over inflationary model parameters. In this case, the table of integrals would 
have to cover the 5D parameter space (a, (3, l\, I2, h) which may seem unrealistic, but our discussion has shown 
otherwise. While the cv/3-triangles must be sampled to fairly high resolution (e.g. on an N = 800 grid), in contrast, 
the CMB bispectrum is intrinsically a slowly varying function of the multipoles (period I ~ 200, like the CMB power 
spectrum). Provided the shape function is relatively featureless, we would only require sparse sampling in multipole 
space (perhaps as few as 33 points on a 2D slice, given the underlying symmetries, even for I ~ 1000). The initial 
tabulation would entail much parallel computation, but subsequently complete bispectra (up to I < 2000) could be 
accurately calculated in minutes. 



Having calculated the bispectrum accurately today we then need to turn to the problem of measurement. The 
bispectrum itself will unfortunately be too small for direct measurement so we will have to rely on estimators. 
Optimal estimators for the bispectrum are based on the statistic, 



where Bi^i^ is the bispectrum predicted from theory and the a; m 's are the measured values. If the normalisation 
factor is, 



then the statistic gives the relative magnitude of the theoretical bispectrum that best fits the data. With the new 
methods described in the previous section we can now calculate the theoretical bispectrum accurately for an almost 
arbitrary shape function. Unfortunately however, we are still be restricted by difficulties in calculating the observed 
bispectrum. One of the biggest challenges is to calculate the Wigner 3j-symbols. These have no simple analytic form 
for the general case and are too many too tabulate for large values of 1. Even for a relatively modest l max — 335 as 
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Higher Derivative Difference 
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Figure 16: Plots of the bispectrum for the higher derivative model (on the left) and the difference with the equilateral case 
(on the right) for I = constant slices using the full radiation transfer function. The top row is 31 = 850, the second 31 = 1850, 
and the third 31 = 3650. The bispectrum has been divided by the analytic result for the local case so its features are clearly 
visible. Lines of constant U are those parallel to their respective edge of the triangle. The extremums of the graphs always 
appear when the U triples contain Us which represent extremums of the plot along the / direction. 

used for analysis of the WMAP 1-year data requires calculating over 80 billion independent 3j-symbols. This makes 
the sum over m impractical to compute even if the aim's were all known. The separable case overcomes this by using 
the Gaunt integral to replacing the 3j-symbol with an integral over three spherical harmonics which are absorbed with 
the aim's into the separable parts allowing for efficient calculation. In the general case this is impossible so we must 
find another way around the problem. 

In a fast method is proposed for calculating the estimator. When the reduced bispectrum can be represented 
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Figure 17: Plots of the bispectmm for the DBI case (on the left) and for the equilateral case (on the right) for / < 1800. Note 
how in the equilateral case the perturbations on the sides are suppressed 



in a separable form, 



Then if we defining, 



we can write the estimator as, 



b hhh = gE { X h Y h Z h + 5 Permutations) . (52) 



X«(n) = ^X«^y im (n), (53) 

lm 



S = M £ / dhX^(h)Y^(h)Z^(h). (54) 
i=l J 

While the general reduced bispectrum has no simple separable form, we note from our earlier plot that it is smooth. 
This means that we should be able to represent it as a sum of smooth basis functions. As we have seen from previous 
plots (fTTl [T2l [TB]) when we divide the bispectrum by the analytic result we get a surface which is 0(1) with slow 
oscillations. As a result we choose our basis as follows, 



b hhh = g ^2 a a/3 7 (X' a (li)Xp(l2)X^(h) + 2 permutations) , (55) 
where X a are shifted Legendre polynomials, 

X a (l) = P a ( 2l 7 lmax ) (56) 



and, 



However, any set of orthogonal polynomials would suffice. Simplifying we have 



h(h + l) + h(l2 + l)+h(h + l) 



7T j hihh = ^ ag^Xg X fi {h)X n (h) (58) 



a/37 
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And a a /3j can be easily determined, 

to _i_iVofl-LiVo , n f d hdl 2 dl 3 ( h(h + ^)h{h + 1)^3 + 1) , , v ,, v , , , r , n 

a Q/37 = (2a + l)(2/3 + 1)(2 7 + 1) / — , ti Kl 3 l 3 X a {li)X^l2)X^h), (59) 

using the orthogonality condition, 



If we define as before. 



and 



dxP a (x)P,(x) = => f^' JL Xa (l)X,(l) = -^-. (60) 

2a + 1 Jr. l max 2a + 1 



X a (n) =J2Xa(l)^Y lm (&), (61) 



Zm 



Then the estimator becomes, 



Where we have defined 



S = — ^ a Q( 3 7 M a/ 3 7 . (63) 

M Q/ 3 7 = - J dn (X' a {p)X' fi {n)X 7 {a) + 2 permutations) . (64) 

This approach has several advantages. Firstly, there is complete separation of theory and measurement. The quantities 
X,X' and M Q/ 3 7 only need to be calculated once per map then M a p-y can simply be stored on disk. For each theory 
we only need to calculate a a( 3 7 then perform the sum of the product of it with the waiting M ai g 7 . This process will 
be quick as rather than summing over all combinations of i, we now only need to do so over the range of a, 0, 7. This 
new approach for the estimator makes it possible to produce quick and accurate bispectrum estimates for general 
models. This decomposition has been tested for l max — 400 for the local model. Using basis functions up to a = 30 
we managed to reconstruct the reduced bispectrum to a 1% accuracy as shown by the plots (fT8l fl"9|) . 



FUTURE WORK 



With the numerical methods we have developed it is now possible to accurately calculate the CMB bispectrum 
today from a general shape function. In the stringent test case there is a problem obtaining accuracy without long 
calculation times because of truncation errors, but arbitrary accuracy can be achieved by altering truncation thresholds 
and improving resolution. In the real calculation using the full radiation transfer function we have demonstrated that 
the systematic error sources are much weaker and very accurate results can be obtained rapidly using our adaptive 
method. We find that the bispectra for the local and equilateral cases are broadly similar in the equal multipole case, 
with the main difference being that the equilateral case is smaller at low I. Where the two cases differ significantly is in 
the slices ^1 + ^2 + ^3 = contant (and for which we find no significant new sources of error). For the calculation with the 
full radiation transfer functions, we find that the results are much the same in the centre, but with the equilateral case 
more heavily suppressed where the / values differ significantly, as expected. This indicates that the equal multipole 
bispectrum provides a good measure of the overall size of the non-Gaussianity, while the unequal I bispectrum provides 
a differentiator between competing models. Obviously, meaningful plots of the CMB bispectrum (like figure H"2"]) are 
not expected to be determined observationally, not least because they will be dominated by noise and convolved with 
experimental beams. Instead we must use an estimator to determine the fit of theory to observation. 

While a parallel implementation of the code runs sufficiently rapidly on a supercomputer for calculations of the 
full bispectrum, there is still plenty of scope for further improvements. Essentially all of the calculation time is spent 
evaluating the one-dimensional integrals with the geometric integral taking twice as long as the transfer integral. 
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Figure 18: Plot comparing the equal I reduced bispectrum and the representation of it in terms of basis functions up to a = 30 

Local Legendre 





Figure 19: Plots of the bispectrum for the local model plotted next to the representation of it in terms of basis functions up to 
a = 30. 



To reduce overall evaluation times we then have two options. We can improve the integration methods for the one- 
dimensional integrals or we can reduce the number of points we need to calculate by improving the adaptive algorithm. 
Improving the cell selection criteria in the adaptive algorithm has largely been a long process of trial and error, and 
no doubt this can be continued with further success. However, improving the one-dimensional integration method is 
perhaps the more exciting option. Here, more rapidly convergent numerical methods are possible, but the hope is 
that for the case of the geometric integral we may be able to switch to an analytic solution should a stable method 
of evaluation be discovered. Such an advance could reduce overall bispectrum calculation times by a factor of 3. 

The decomposition of the bispectrum into Legendre polynomials has also been parallelised but is still hampered but 
the requirement to use a large number of polynomials. This arises from having to perform the integral over the cube 
in 1 space with sides from to l ma x while the bispectrum is only defined inside the tetrahedron created by the triangle 
condition on the Vs. When we use the code to calculate points outside this tetrahedron we introduce discontinuities 
which require a large number of polynomials to fit accurately. We are currently looking to edit the code to calculate 
these points so that the bispectrum remains smooth as we cross the boundary. 

The present code has been developed in modular form to link and work together with existing CMB line-of-sight 
codes for the angular power spectrum (like CMBFAST, CAMB or CMB2000); it extracts the appropriate radiation 
transfer functions from these for a given cosmology. In due course, it is our intention to improve the efficiency, 
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portability and parallelism (as well as the documentation) of this code and to make it publicly available. 

The bispectrum remains a highly significant tool for testing for non-Gaussianity in the CMB. Hopefully, future 
experiments like the Planck satellite will either detect or provide much tighter constraints on this non-Gaussianity, 
opening a new window on the early universe. With this code we now have the capability to make accurate predictions 
for the contributions to the CMB bispectrum which are of primordial origin. For example, we plan to use quantitative 
bispectra from multifield inflation models [|| and project these forward to make falsifiable CMB predictions. Clearly 
parallel efforts must also continue in order to pin down and distinguish alternative sources of CMB non-Gaussianity 
from later times, such as second-order effects in the radiation transfer functions and nonlinear astrophysical effects. 
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THREE BESSEL FUNCTION INTEGRALS 



Higher-order angular correlation functions generically involve integrals over products of Bessel functions. If we 
allow a general bispectrum but remain in the large angle approximation then both the ID integrals in (|28p . I T (a,(3) 
and I G (a, (3), are geometric and of this form. Assuming n = 0, 



1 



dk 



1 ( a ,P) = TvT 3h(ak)ji 2 (bk)ji 3 (ck)—, 



3 3 ./ •' " J ' ' k 

where Ar] has been absorbed via a rescaling of k. This leaves us trying to solve integrals of the form, 



Ih,h,l 3 (a,b, c,n) 



3h ( ax ) ih ( bx ) 3h ( cx ) xtldx 



(65) 



(66) 



where n = 2, —1. 

We already have a several solutions for n = 2 given in refs. 15, 2, 17 1. For the general case a solution is proposed 
in ref. [2~i| but, unfortunately, their analysis is incomplete. Their method revolves around replacing the spherical 
Bcsscl functions as below 
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to rewrite the integral as, 
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The second line of the integrand can be expanded to, 
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So we are trying to complete the integral [68] for eight values of I, 



dx. 



(70) 



This integral is singular but the sum is not. The authors calculate, 

T dx, 



(x + e) p 

and assert that the infinite terms cancel explicitly in the sum as e — > 0, which is incorrect. If we instead calculate 



(71) 
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dx. 



which has solution, 



(il) 



p-i 



p-2 



(P-1) 

where E\ is the En- Function with series expansion, 
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Eri-iel) = - 7 - ln(e) - In \l\ + -sgn(l) + 0(e), 



(72) 



(73) 



(74) 



and 7 is the Euler constant. This is identical to that in ref. 24j result except for the e lei factor in front of the sum. 
The claim that the infinite part cancelled led the authors to drop the sum altogether. Here the singular parts of the 
sum over s cancel in the triple sum over ra„ and as we will set e to zero, we are then only left with in the constant 
part of the series. 
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Where H a = 2^s=i ^ s the harmonic number. 

In the triple sum over m, the parts proportional to 7 and ln(e) in E\ cancel exactly. The part of the integral that 
contributes to the final result is then, 
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We can then substitute this into the expression for the integral over the spherical Bessel functions to eventually obtain, 
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where, 

(a + b + c)™i+™'+ m 3-n+2 (gmi+m2+m3 _ n+2 _ log | Q + 6 + c |) 

(2a) m i+ 1 (26)™2+i(2 c ) m 3+i [ ' S > 

(a + 6 + c) mi + m2 + m3 -™+ 2 sign(a+b+c) 

(2a) m i+ 1 (26) m 2+i(2c) m 3+i ' ^ 79 ) 

These analytic solutions have unfortunately not proved particularly useful in the evaluation of the bispectrum. The 
solutions that involve large series like above are unstable for large I's and also in the corners of the a/3-triangles, where 
one of the k^s tends to zero. The difficulty with these series for a straightforward numerical implementation is that 
they generically involve the exact cancellation of many very large terms. 
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